Installation and requirements

serosolver may be installed from github using the devtools package. There are a number of additional packages that we need for this analysis.

# Required to run serosolver
#devtools::install_github("seroanalytics/serosolver")
library(serosolver)
library(plyr)
library(data.table)

## Required for this analysis
library(reshape2)
library(foreach)
library(doParallel)
library(bayesplot)
library(coda)
library(ggplot2)
library(viridis)

## doRNG handles seeding in the foreach loop
## Note that this vignette was generated on a Linux machine,
## and the setup for parallelisation is different on a Windows machine
set.seed(1234)
cl <- makeCluster(5)
registerDoParallel(cl)

## Choose wd to run anaylsis
setwd("~/Documents/serosolver/vignettes")

Assumptions

In this analysis, serological samples were taken between 2009 and 2011 and therefore all time variables are relative to this time period. We are interested in inferring infections and attack rates at a quarterly resolution, and therefore set resolution to 4. Our primary outcome of interest is to infer unbiased attack rates, and we therefore use the version of the code with a Beta prior on per-time attack rates, prior_version=2. We set these parameters at the start of the analysis. Additionally, we assume that the samples are tested against the same virus, therefore, we set all antigenic coordinates (more information below) to those of the first time point.

filename <- "case_study_1"
resolution <- 4 ## set to 4 for quarterly resolution
sample_years <- 2009:2012

serosolver::describe_proposals()
#> [1] "Which version to use in run_MCMC? The following text describes the proposal step for updating infection histories."
#> [1] "Version 1: Beta prior on per time attack rates. Explicit FOI on each epoch using probability of infection term. Proposal performs N `flip` proposals at random locations in an individual's infection history, switching 1->0 or 0->1. Otherwise, swaps the contents of two random locations"
#> [1] "Version 2: Beta prior on per time attack rates. Gibbs sampling of infection histories as in Indian Buffet Process papers, integrating out each probability of infection term."
#> [1] "Version 3: Beta prior on probability of infection for an individual, assuming independence between individuals. Samples from a beta binomial with alpha and beta specified by the par_tab input. Proposes nInfs moves at a time for add/remove, or when swapping, swaps locations up to moveSize time steps away"
#> [1] "Version 4: Beta prior on probability of any infection. Gibbs sampling of infection histories using total number of infections across all times and all individuals as the prior"
prior_version <- 2

Preparing the data

The data used in this analysis are haemagglutination inhibition (HI) titres against A/H1N1pdm09 that began circulating in 2009. The raw data have been pre-processed to both convert them into a form usable for serosolver and to separate the data into vaccinated and unvaccinated data sets. Details about the pre-processing can be found in the scripts folder of the serosolver package.

## Read in titre data
 inputDat <- read.csv(file = "case_study_1/data/HKdata_h1n1_unvac.csv",header = TRUE)
# inputDat <- read.csv(file = "case_study_1/data/HKdata_h1n1_vac.csv",header = TRUE)

indivs <- unique(inputDat$individual) #all individuals

# Subset data for indivs
titre_dat <- inputDat[inputDat$individual %in% indivs,c("individual","virus","titre","samples","DOB")]
titre_dat$individual <- match(titre_dat$individual, indivs)

titre_dat <- unique(titre_dat)
titre_dat <- plyr::ddply(titre_dat,.(individual,virus,samples),function(x) cbind(x,"run"=1:nrow(x),"group"=1))
print(head(titre_dat))
#>   individual virus titre samples  DOB run group
#> 1          1  8037     0    8039 8036   1     1
#> 2          1  8037     0    8040 8036   1     1
#> 3          1  8037     7    8044 8036   1     1
#> 4          1  8037     7    8047 8036   1     1
#> 5          2  8037     0    8039 8036   1     1
#> 6          2  8037     5    8041 8036   1     1

NOTE: vaccinated and unvaccinated data must be run separately.

Given that this analysis uses titres from a single virus it is necessary to define an antigenic map with matching antigenic distance for each time point. We use coordinates based on the antigenic map created by Fonville et al. Generating the antigenic map involves fitting a smoothing spline through provided coordinates to give a representative virus for each time point (in this case, each quarter) that an individual could be infected. This process also imputs antigenic coordinates for time points that we do not have a measured virus. We then overwrite the coordinates at every time point with those of the first time point.

## Read in raw coordinates
antigenic_coords <- read.csv("case_study_1/data/fonville_map_approx.csv",stringsAsFactors=FALSE)

## Convert to form expected by serosolver
antigenic_map <- generate_antigenic_map(antigenic_coords, resolution)

## Restrict entries to years of interest. Entries in antigenic_map determine
## the times that individual can be infected ie. the dimensions of the infection
## history matrix.
antigenic_map <- antigenic_map[antigenic_map$inf_years>=(sample_years[1]*resolution+1) & antigenic_map$inf_years<=sample_years[4]*resolution,]

## Change all coordinates to the same as the first pair of coordinates (because we assume the virus is the same at every time point)
antigenic_map[2:nrow(antigenic_map),c('x_coord','y_coord')] <- antigenic_map[1,c('x_coord','y_coord')]
print(head(antigenic_map))
#>      x_coord  y_coord inf_years
#> 166 37.69244 8.233511      8037
#> 167 37.69244 8.233511      8038
#> 168 37.69244 8.233511      8039
#> 169 37.69244 8.233511      8040
#> 170 37.69244 8.233511      8041
#> 171 37.69244 8.233511      8042

strain_isolation_times <- unique(antigenic_map$inf_years)

NOTE: generate_antigenic_map expects the provided file fonville_map_approx.csv. Users should refer to generate_antigenic_map_flexible for more generic antigenic map generation.

Finally, we must specify the par_tab data frame, which controls which parameters are included in the model, which are fixed, and their uniform prior ranges. Given that we are integrating out the probability of infection terms under prior version 2, we must remove these parameters from par_tab. Furthermore, given that we are interested in short-term dynamics with relatively sparse data, we remove parameters relating to the long-term antibody kinetics phase to avoid identifiability issues. We set alpha and beta to 1.

par_tab <- read.csv("case_study_1/data/parTab_base.csv",stringsAsFactors=FALSE)

## Set parameters for beta and alpha to 1
par_tab[par_tab$names %in% c("alpha","beta"),"values"] <- c(1,1)
## Maximum recordable log titre in these data is 9
par_tab[par_tab$names == "MAX_TITRE","values"] <- 9

## Remove phi parameters, as these are integrated out under prior version 2
par_tab <- par_tab[par_tab$names != "phi",]

## Fix all long term parameters to 0
par_tab[par_tab$names %in% c("mu","tau","sigma1","sigma2"),"fixed"] <- 1 #mu, tau, sigma1, and sigma2 are fixed
par_tab[par_tab$names %in% c("mu","tau","sigma1","sigma2"),"values"] <- 0 # set these values to 0

Summary

• Choose resolution, attack rate priors and reference time for “the present”
• Read in titre data and convert into a useable form for serosolver
• Generate antigenic map for the exposing virus
• Generate parameter control table for MCMC

Extract attack rates from fits

We simulate infection histories and antibody titre data based on the “real” parameters inferred from fitting the model above. First, we extract the maximum posterior probability antibody kinetics parameters and attack rates.

## Read in MCMC chains from fitting
all_chains <- load_mcmc_chains(location="case_study_1/chains_main/",thin=10,burnin=10000,
                               par_tab=par_tab,unfixed=FALSE,convert_mcmc=FALSE)
#> Chains detected:     5[[1]]
#> [1] 501
#> 
#> [[2]]
#> [1] 501
#> 
#> [[3]]
#> [1] 501
#> 
#> [[4]]
#> [1] 501
#> 
#> [[5]]
#> [1] 501
#> 
#> Chains detected: 
#> case_study_1/chains_main//case_study_1_1_infection_histories.csv
#> case_study_1/chains_main//case_study_1_2_infection_histories.csv
#> case_study_1/chains_main//case_study_1_3_infection_histories.csv
#> case_study_1/chains_main//case_study_1_4_infection_histories.csv
#> case_study_1/chains_main//case_study_1_5_infection_histories.csv
#> [[1]]
#> [1] 77778
#> 
#> [[2]]
#> [1] 77622
#> 
#> [[3]]
#> [1] 77613
#> 
#> [[4]]
#> [1] 77432
#> 
#> [[5]]
#> [1] 77549

## Find samples that were in both theta and inf hist chains
chain <- all_chains$theta_chain
inf_chain <- all_chains$inf_chain
intersect_samps <- intersect(unique(inf_chain$sampno), unique(chain$sampno))
chain <- chain[chain$sampno %in% intersect_samps,]

## Find the parameter values that gave the highest posterior probability
which_mle <- chain[which.max(chain$lnlike),c("sampno","chain_no")]
mle_theta_pars <- chain[chain$sampno == which_mle$sampno & chain$chain_no == which_mle$chain_no,]

## Store total infections to compare later
mle_total_infs <- mle_theta_pars[,"total_infections"]
mle_theta_pars <- mle_theta_pars[,par_tab$names]
mle_inf_hist <- inf_chain[inf_chain$sampno == which_mle$sampno & inf_chain$chain_no == which_mle$chain_no,]

## Generate full infection history matrix using provided function
mle_inf_hist <- expand_summary_inf_chain(mle_inf_hist[,c("sampno","j","i","x")])
## Find number of infections per year from this infection history
no_infs <- colSums(mle_inf_hist[,3:ncol(mle_inf_hist)])

## If missing time points in simulated attack rates
if(length(no_infs) < length(strain_isolation_times)){
  diff_lengths <- length(strain_isolation_times) - length(no_infs)
  no_infs <- c(no_infs, rep(0, diff_lengths))
}

## Find attack rate per year
n_alive <- get_n_alive(titre_dat, strain_isolation_times)
attack_rates <- no_infs/n_alive

Functions are provided to simulate antibody titre data under a given serosurvey design. The antibody kinetics parameters and attack rates estimated above are used to simulate titres from the model. The simulate_data function is well documented, and users should refer to the help file to customise the simulated serosurvey design.

set.seed(1234)

sim_par_tab <- par_tab
sim_par_tab$values <- as.numeric(mle_theta_pars)
sim_par_tab[sim_par_tab$names %in% c("alpha","beta"),"values"] <- c(1/3,1/3)
sim_par_tab[sim_par_tab$names == "wane","values"] <- 1
sim_par_tab[sim_par_tab$names == "wane","values"] <- par_tab[par_tab$names == "wane","values"]/resolution
sim_par_tab[sim_par_tab$names %in% c("mu","tau","sigma1","sigma2"),"fixed"] <- 1 
sim_par_tab[sim_par_tab$names %in% c("mu","tau","sigma1","sigma2"),"values"] <- 0 
sim_par_tab[sim_par_tab$names == "MAX_TITRE","values"] <- 9 

sampling_times <- seq(2009*resolution + 1, 2012*resolution, by=1)

age_min <- 6*resolution
age_max <- 6*resolution
n_indiv <- length(unique(titre_dat$individual))

dat <- simulate_data(par_tab = sim_par_tab, 
                     n_indiv = n_indiv,
                     buckets = resolution,
                     strain_isolation_times = strain_isolation_times,
                     sampling_times = sampling_times, 
                     nsamps = 4, 
                     antigenic_map = antigenic_map, 
                     age_min = age_min,
                     age_max = age_max,
                     attack_rates=attack_rates,
                     repeats = 1)



## Inspect simulated antibody titre data and infection histories
sim_titre_dat <- dat[["data"]]
sim_infection_histories <- dat[["infection_histories"]]

## Store total infections to compare later
actual_total_infections <- sum(sim_infection_histories)

plot_data(sim_titre_dat, sim_infection_histories, strain_isolation_times,n_samps = 5)

## Use titres only against same viruses tested in real data
viruses <- unique(titre_dat$virus)
sim_titre_dat <- sim_titre_dat[sim_titre_dat$virus %in% viruses, ]
sim_ages <- dat[["ages"]]
sim_titre_dat <- merge(sim_titre_dat, sim_ages)
sim_ar <- dat[["attack_rates"]]

Simulation fitting

Once these simulated data have been generated, the work flow becomes exactly the same as with the real data above.

filename <- "case_study_1_sim"

## Distinct filename for each chain
no_chains <- 5
filenames <- paste0(filename, "_",1:no_chains)

## Create the posterior solving function that will be used in the MCMC framework 
f <- create_posterior_func_fast(par_tab=sim_par_tab,
                                titre_dat=sim_titre_dat,
                                antigenic_map=antigenic_map,
                                version=prior_version) # function in posteriors.R

## Generate results in parallel
res <- foreach(x = filenames) %dopar% {
  ## Not all random starting conditions return finite likelihood, so for each chain generate random
  ## conditions until we get one with a finite likelihood
  start_prob <- -Inf
  while(!is.finite(start_prob)){
    start_tab <- sim_par_tab
    for(i in 1:nrow(start_tab)){
      if(start_tab[i,"fixed"] == 0){
        start_tab[i,"values"] <- runif(1,start_tab[i,"lower_start"], 
                                      start_tab[i,"upper_start"])
      }
    }
    ## Generate starting infection history
    start_inf <- setup_infection_histories_new_2(sim_titre_dat, strain_isolation_times, space=3,titre_cutoff=4)
    start_prob <- sum(f(start_tab$values, start_inf))
  }
  
  res <- run_MCMC(par_tab = start_tab, 
                  titre_dat = sim_titre_dat,
                  antigenic_map = antigenic_map,
                  start_inf_hist = start_inf, 
                  mcmc_pars = c("iterations"=500000,"adaptive_period"=100000,"thin"=10,
                                "thin_hist"=1000,"save_block"=1000,
                                "hist_switch_prob"=0.8, "hist_sample_prob"=1),
                  filename = paste0("chains_sim/",x), 
                  CREATE_POSTERIOR_FUNC = create_posterior_func_fast, 
                  version = prior_version,
                  fast_version=TRUE)
}

Simulation analysis

MCMC chains should be checked for convergence under the usual diagnostics. We also compare the inferred posterior distributions to the known true parameter values. We see that convergence and between-chain agreement is good and that the model recovers reasonably unbiased estimates for some parameters. However, under this sampling strategy the model slightly underestimates the amount of long term antibody boosting elicited by a single infection and overestimates the total number of infections. This is driven by the contribution of the attack rate prior relative to the contribution of the likelihood (the data). Increasing the number of measured titres (for example, measure titres against 40 viruses rather than 9) or using a more informative attack rate prior would help reduce this bias.

## Read in the MCMC chains
## Note that `thin` here is in addition to any thinning done during the fitting
sim_all_chains <- load_mcmc_chains(location="case_study_1/chains_sim",thin=1,burnin=10000,
                                   par_tab=par_tab,unfixed=FALSE,convert_mcmc=TRUE)
#> Chains detected:     5[[1]]
#> [1] 45001
#> 
#> [[2]]
#> [1] 45001
#> 
#> [[3]]
#> [1] 45001
#> 
#> [[4]]
#> [1] 45001
#> 
#> [[5]]
#> [1] 45001
#> 
#> Chains detected: 
#> case_study_1/chains_sim/case_study_1_sim_1_infection_histories.csv
#> case_study_1/chains_sim/case_study_1_sim_2_infection_histories.csv
#> case_study_1/chains_sim/case_study_1_sim_3_infection_histories.csv
#> case_study_1/chains_sim/case_study_1_sim_4_infection_histories.csv
#> case_study_1/chains_sim/case_study_1_sim_5_infection_histories.csv
#> [[1]]
#> [1] 597453
#> 
#> [[2]]
#> [1] 596783
#> 
#> [[3]]
#> [1] 596893
#> 
#> [[4]]
#> [1] 597311
#> 
#> [[5]]
#> [1] 596685
theta_chain <- sim_all_chains$theta_chain
## Get the MCMC chains as a list
list_chains <- sim_all_chains$theta_list_chains
## Look at diagnostics for the free parameters
list_chains1 <- lapply(list_chains, function(x) x[,c("mu_short","wane","error",
                                                     "total_infections",
                                                     "lnlike","prior_prob")])

## Gelman-Rubin diagnostics and effective sample size
print(gelman.diag(as.mcmc.list(list_chains1)))
#> Potential scale reduction factors:
#> 
#>                  Point est. Upper C.I.
#> mu_short               1.00       1.01
#> wane                   1.00       1.01
#> error                  1.00       1.01
#> total_infections       1.00       1.00
#> lnlike                 1.03       1.08
#> prior_prob             1.03       1.07
#> 
#> Multivariate psrf
#> 
#> 1.04
print(effectiveSize(as.mcmc.list(list_chains1)))
#>         mu_short             wane            error total_infections 
#>        1875.7357        1598.0157        3427.4839        4027.1533 
#>           lnlike       prior_prob 
#>         418.7965         489.4283

melted_theta_chain <- reshape2::melt(as.data.frame(theta_chain), id.vars=c("sampno","chain_no"))
estimated_pars <- c(sim_par_tab[sim_par_tab$fixed == 0,"names"],"total_infections")
melted_theta_chain <- melted_theta_chain[melted_theta_chain$variable %in% estimated_pars,]
colnames(melted_theta_chain)[3] <- "names"

add_row <- data.frame("total_infections",actual_total_infections,0,0.1,0,10000,0,0,1)
colnames(add_row) <- colnames(sim_par_tab)
sim_par_tab1 <- rbind(sim_par_tab, add_row)

ggplot(melted_theta_chain) + 
  geom_density(aes(x=value,fill=as.factor(chain_no)),alpha=0.5) +
  geom_vline(data=sim_par_tab1[sim_par_tab1$fixed == 0,],aes(xintercept=values),linetype="dashed") +
  facet_wrap(~names,scales="free") + 
  theme_classic() +
  theme(legend.position="bottom")

Recovery of known attack rates is also reasonably accurate, though the constraint of the posterior distibution is quite low for many years where identifiability is poor. Again, more titre data or more individuals would improve inferential power. One particularly reassuring plot is the comparison of known individual cumulative infection histories (the cumulative sum of infections over time for an individual) against the estimated posterior distribution of cumulative infection histories. We see that the 95% credible intervals capture the true cumulative infection histories in almost all cases.

## Extract infection history chain
inf_chain <- sim_all_chains$inf_chain

## Look at inferred attack rates
p_ar <- plot_attack_rates(inf_chain, sim_titre_dat, strain_isolation_times, pad_chain=FALSE,
                            plot_den = TRUE,prior_pars=list(prior_version=prior_version, 
                                                          alpha=par_tab[par_tab$names=="alpha","values"],
                                                          beta=par_tab[par_tab$names=="beta","values"]))  +
  geom_point(data=sim_ar,aes(x=year,y=AR),col="purple")
print(p_ar)

## Calculate convergence diagnostics and summary statistics on infection histories
## Important to scale all infection estimates by number alive from titre_dat
sim_n_alive <- get_n_alive(sim_titre_dat, strain_isolation_times)

## This function generates a number of MCMC outputs
ps_infhist <- plot_posteriors_infhist(inf_chain=inf_chain, 
                                      years=strain_isolation_times, 
                                      n_alive=sim_n_alive)

## Check for agreement between inferred cumulative infection histories 
## for some individuals
p_indiv_inf_hists <- generate_cumulative_inf_plots(inf_chain,indivs=c(12,21,39,41,44,66,71,76,79),pad_chain=FALSE,
                                                   real_inf_hist=sim_infection_histories,
                                                  strain_isolation_times = strain_isolation_times,
                                                  number_col=3)
print(p_indiv_inf_hists[[1]])

## Posterior probability that infections occured at given times per individual
print(p_indiv_inf_hists[[2]])